#include "defines.h"


int test(DATA_TYPE * B, int NX, int NY, int T, int Bx, int By, int tb) {

	
	if(Bx<(2*XSLOPE+1) || By<(2*YSLOPE+1) || Bx>NX || By>NY || tb>min(((Bx-1)/2)/XSLOPE,((By-1)/2)/YSLOPE)){
		return 0;
	}
	
	DATA_TYPE (*A)[NX+2*XSTART][NY+2*YSTART] = (DATA_TYPE (*)[NX+2*XSTART][NY+2*YSTART]) B;


	int bx = Bx-2*(tb*XSLOPE);
	int by = By-2*(tb*YSLOPE);

	int ix = Bx+bx;
	int iy = By+by; 

	int xnb0  = ceild(NX,ix);
	int ynb0  = ceild(NY,iy);
	int xnb11 = ceild(NX+(Bx-bx)/2,ix); // the start x dimension of B11 is bx,
	int ynb12 = ceild(NY+(By-by)/2,iy); // for periodic boundary stencils xnb11 and ynb12 must be larger than 2.
	int ynb11 = ynb0;
	int xnb12 = xnb0;
	int xnb2  = xnb11;
	int ynb2  = ynb12;

	int nb1[2]   = {xnb11 * ynb11, xnb12 * ynb12};
	int nb02[2]  = {xnb0 * ynb0, xnb2 * ynb2};
	int xnb1[2]  = {xnb11, xnb12};
	int xnb02[2] = {xnb0, xnb2};

	int xleft02[2]   = {XSTART + bx, XSTART - (Bx-bx)/2};
	int ybottom02[2] = {YSTART + by, YSTART - (By-by)/2};
	int xleft11[2]   = {XSTART,		 XSTART + ix/2};
	int ybottom11[2] = {YSTART + by, YSTART - (By-by)/2};
	int xleft12[2]   = {XSTART + bx, XSTART - (Bx-bx)/2};
	int ybottom12[2] = {YSTART,		 YSTART + iy/2};

	int level = 0;
	int tt,n;
	int x, y, t;
	register int ymin, ymax;
	int xmin,xmax;



	for(tt =- tb; tt < T; tt += tb){
#pragma omp parallel for schedule(dynamic)  private(xmin,xmax,ymin,ymax,t,x,y)
		for(n = 0; n < nb02[level]; n++){

			for(t = max(tt,0); t < min(tt + 2*tb, T); t++) { 

				xmin = max(   XSTART,   xleft02[level] + (n%xnb02[level]) * ix      + myabs(t+1,tt+tb) * XSLOPE);
				xmax = min(NX+XSTART,   xleft02[level] + (n%xnb02[level]) * ix + Bx - myabs(t+1,tt+tb) * XSLOPE);
				ymin = max(   YSTART, ybottom02[level] + (n/xnb02[level]) * iy      + myabs(t+1,tt+tb) * YSLOPE);
				ymax = min(NY+YSTART, ybottom02[level] + (n/xnb02[level]) * iy + By - myabs(t+1,tt+tb) * YSLOPE);

				for(x = xmin; x < xmax; x++) {
#pragma ivdep
#pragma vector always
					for(y = ymin; y < ymax; y++) {
						kernel(A);
					}
				}
			}
		}


#pragma omp parallel for schedule(dynamic)  private(xmin,xmax,ymin,ymax,t,x,y)
		for(n = 0; n < nb1[0] + nb1[1]; n++) {

			for(t = tt + tb; t < min(tt + 2*tb, T); t++) {
				if(n < nb1[level]) {
					xmin = max(     XSTART,   xleft11[level] + (n%xnb1[level]) * ix       - (t+1-tt-tb) * XSLOPE);
					xmax = min(NX + XSTART,   xleft11[level] + (n%xnb1[level]) * ix  + bx + (t+1-tt-tb) * XSLOPE);
					ymin = max(     YSTART, ybottom11[level] + (n/xnb1[level]) * iy       + (t+1-tt-tb) * YSLOPE);
					ymax = min(NY + YSTART, ybottom11[level] + (n/xnb1[level]) * iy  + By - (t+1-tt-tb) * YSLOPE);
				}
				else {

					xmin = max(     XSTART,   xleft12[level] + ((n-nb1[level])%xnb1[1-level]) * ix      + (t+1-tt-tb) * XSLOPE);
					xmax = min(NX + XSTART,   xleft12[level] + ((n-nb1[level])%xnb1[1-level]) * ix + Bx - (t+1-tt-tb) * XSLOPE);
					ymin = max(     YSTART, ybottom12[level] + ((n-nb1[level])/xnb1[1-level]) * iy      - (t+1-tt-tb) * YSLOPE);
					ymax = min(NY + YSTART, ybottom12[level] + ((n-nb1[level])/xnb1[1-level]) * iy + by + (t+1-tt-tb) * YSLOPE);

				}
				for(x = xmin; x < xmax; x++) {
#pragma ivdep
#pragma vector always
					for(y = ymin; y < ymax; y++) {
						kernel(A);
					}
				}
			}
		}
		level = 1- level;
	}
}
